Background

Synthetic data assimilation experiments known as perfect model assimilation are a convenient way to “know” the truth of what is being assimilated. This can help evaluate the data assimilation system and the limits of skill improvements through DA.

Setup and Visualize

Load the rwrfhydro package.

library("rwrfhydro")
## To check rwrfhydro updates run: CheckForUpdates()
# set your own path
rlFile <- '~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03//Route_Link.nc' 
# currently there is no "gages" variable
ncdump(rlFile)  
## File: ~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03//Route_Link.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (1):
##     linkDim = 403 ; 
## variables (24):
##     int link(linkDim)    ; 
##         link:long_name = "Link ID" ;
##     int from(linkDim)    ; 
##         from:long_name = "From Link ID" ;
##         from:coordinates = "time lat lon" ;
##     int to(linkDim)    ; 
##         to:long_name = "To Link ID" ;
##         to:coordinates = "time lat lon" ;
##     float lon(linkDim)    ; 
##         lon:long_name = "longitude of the start node" ;
##         lon:units = "degrees_east" ;
##         lon:standard_name = "longitude" ;
##     float lat(linkDim)    ; 
##         lat:long_name = "latitude of the start node" ;
##         lat:units = "degrees_north" ;
##         lat:standard_name = "latitude" ;
##     float alt(linkDim)    ; 
##         alt:long_name = "Elevation in meters at start node" ;
##         alt:standard_name = "height" ;
##         alt:units = "m" ;
##         alt:positive = "up" ;
##         alt:axis = "Z" ;
##     int type(linkDim)    ; 
##         type:long_name = "Link type" ;
##         type:coordinates = "time lat lon" ;
##     int order(linkDim)    ; 
##         order:long_name = "Stream order (Strahler)" ;
##         order:coordinates = "time lat lon" ;
##     float Qi(linkDim)    ; 
##         Qi:long_name = "Initial flow in link (CMS)" ;
##         Qi:coordinates = "time lat lon" ;
##     float MusK(linkDim)    ; 
##         MusK:long_name = "Muskingum routing time (s)" ;
##         MusK:coordinates = "time lat lon" ;
##     float MusX(linkDim)    ; 
##         MusX:long_name = "Muskingum weighting coefficient" ;
##         MusX:coordinates = "time lat lon" ;
##     float Length(linkDim)    ; 
##         Length:long_name = "Stream length (m)" ;
##         Length:coordinates = "time lat lon" ;
##     float n(linkDim)    ; 
##         n:long_name = "Manning's roughness" ;
##         n:coordinates = "time lat lon" ;
##     float So(linkDim)    ; 
##         So:long_name = "Slope (%; drop/length)" ;
##         So:coordinates = "time lat lon" ;
##     float ChSlp(linkDim)    ; 
##         ChSlp:long_name = "Channel side slope (%; drop/length)" ;
##         ChSlp:coordinates = "time lat lon" ;
##     float BtmWdth(linkDim)    ; 
##         BtmWdth:long_name = "Bottom width of channel" ;
##         BtmWdth:coordinates = "time lat lon" ;
##     float LkHZArea(linkDim)    ; 
##         LkHZArea:long_name = "Maximum lake area (sq. m)" ;
##         LkHZArea:coordinates = "time lat lon" ;
##     float LkMxH(linkDim)    ; 
##         LkMxH:long_name = "Maximum lake elevation (m)" ;
##         LkMxH:coordinates = "time lat lon" ;
##     float WeirC(linkDim)    ; 
##         WeirC:long_name = "Weir coefficient" ;
##         WeirC:coordinates = "time lat lon" ;
##     float WeirL(linkDim)    ; 
##         WeirL:long_name = "Weir length (m)" ;
##         WeirL:coordinates = "time lat lon" ;
##     float OrificeC(linkDim)    ; 
##         OrificeC:long_name = "Orifice coefficient" ;
##         OrificeC:coordinates = "time lat lon" ;
##     float OrificeA(linkDim)    ; 
##         OrificeA:long_name = "Orifice croess-sectional area (sq. m)" ;
##         OrificeA:coordinates = "time lat lon" ;
##     float OrificeE(linkDim)    ; 
##         OrificeE:long_name = "Orifice elevation (m)" ;
##         OrificeE:coordinates = "time lat lon" ;
##     float time()    ; 
##         time:standard_name = "time" ;
##         time:long_name = "time of measurement" ;
##         time:units = "days since 2000-01-01 00:00:00" ;
## 
## // global attributes (2):
##     :featureType = "point"
##     :history = "Created Thu Sep 03 09:30:05 2015"
# This function returns a function which plots the data
PlotRouteLink <- VisualizeRouteLink(rlFile)
PlotRouteLink() # basic plot
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.940325,-105.430153&zoom=8&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

args(PlotRouteLink) # what are the options? unfortunately I havent really documented these.
## function (location = c(lon = mean(range(rl$lon)), lat = mean(range(rl$lat))), 
##     indices = FALSE, comIds = FALSE, zoom = if (length(gageZoom)) 13 else 8, 
##     source = "google", maptype = "terrain", padPlot = 0.01, linkColor = "blue", 
##     textColor = "darkred", gageColor = "cyan", gageShape = 9, 
##     gageZoom = NULL, plotPath = NULL, plotType = "pdf", doPlot = TRUE) 
## NULL
PlotRouteLink(zoom=10, comIds=TRUE) # links on the edge of the domain arent fully plotted... because there's no "to" lat/lon.
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.940325,-105.430153&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

Create the gages variable

Define the synthetic gages to be used for the domain. In this example put gages on ALL the links.

rlLink <- ncdump(rlFile, 'link', quiet=TRUE)
gageIds <- paste0('g',formatC(rlLink,width=3,flag='0'))
newCopyId <- 'gagesAllLinks'
newFile <- AddRouteLinkGage(rlFile, gageIds, rlLink, new=newCopyId, overwrite=TRUE)
print(newFile)
## [1] "/Users/jamesmcc/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc"
ncdump(newFile)
## File: /Users/jamesmcc/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (2):
##     linkDim = 403 ; 
##     IDLength = 15 ; 
## variables (26):
##     integer IDLength(IDLength) ; 
##         IDLength:long_name = "IDLength"
##     int link(linkDim)    ; 
##         link:long_name = "Link ID" ;
##     int from(linkDim)    ; 
##         from:long_name = "From Link ID" ;
##         from:coordinates = "time lat lon" ;
##     int to(linkDim)    ; 
##         to:long_name = "To Link ID" ;
##         to:coordinates = "time lat lon" ;
##     float lon(linkDim)    ; 
##         lon:long_name = "longitude of the start node" ;
##         lon:units = "degrees_east" ;
##         lon:standard_name = "longitude" ;
##     float lat(linkDim)    ; 
##         lat:long_name = "latitude of the start node" ;
##         lat:units = "degrees_north" ;
##         lat:standard_name = "latitude" ;
##     float alt(linkDim)    ; 
##         alt:long_name = "Elevation in meters at start node" ;
##         alt:standard_name = "height" ;
##         alt:units = "m" ;
##         alt:positive = "up" ;
##         alt:axis = "Z" ;
##     int type(linkDim)    ; 
##         type:long_name = "Link type" ;
##         type:coordinates = "time lat lon" ;
##     int order(linkDim)    ; 
##         order:long_name = "Stream order (Strahler)" ;
##         order:coordinates = "time lat lon" ;
##     float Qi(linkDim)    ; 
##         Qi:long_name = "Initial flow in link (CMS)" ;
##         Qi:coordinates = "time lat lon" ;
##     float MusK(linkDim)    ; 
##         MusK:long_name = "Muskingum routing time (s)" ;
##         MusK:coordinates = "time lat lon" ;
##     float MusX(linkDim)    ; 
##         MusX:long_name = "Muskingum weighting coefficient" ;
##         MusX:coordinates = "time lat lon" ;
##     float Length(linkDim)    ; 
##         Length:long_name = "Stream length (m)" ;
##         Length:coordinates = "time lat lon" ;
##     float n(linkDim)    ; 
##         n:long_name = "Manning's roughness" ;
##         n:coordinates = "time lat lon" ;
##     float So(linkDim)    ; 
##         So:long_name = "Slope (%; drop/length)" ;
##         So:coordinates = "time lat lon" ;
##     float ChSlp(linkDim)    ; 
##         ChSlp:long_name = "Channel side slope (%; drop/length)" ;
##         ChSlp:coordinates = "time lat lon" ;
##     float BtmWdth(linkDim)    ; 
##         BtmWdth:long_name = "Bottom width of channel" ;
##         BtmWdth:coordinates = "time lat lon" ;
##     float LkHZArea(linkDim)    ; 
##         LkHZArea:long_name = "Maximum lake area (sq. m)" ;
##         LkHZArea:coordinates = "time lat lon" ;
##     float LkMxH(linkDim)    ; 
##         LkMxH:long_name = "Maximum lake elevation (m)" ;
##         LkMxH:coordinates = "time lat lon" ;
##     float WeirC(linkDim)    ; 
##         WeirC:long_name = "Weir coefficient" ;
##         WeirC:coordinates = "time lat lon" ;
##     float WeirL(linkDim)    ; 
##         WeirL:long_name = "Weir length (m)" ;
##         WeirL:coordinates = "time lat lon" ;
##     float OrificeC(linkDim)    ; 
##         OrificeC:long_name = "Orifice coefficient" ;
##         OrificeC:coordinates = "time lat lon" ;
##     float OrificeA(linkDim)    ; 
##         OrificeA:long_name = "Orifice croess-sectional area (sq. m)" ;
##         OrificeA:coordinates = "time lat lon" ;
##     float OrificeE(linkDim)    ; 
##         OrificeE:long_name = "Orifice elevation (m)" ;
##         OrificeE:coordinates = "time lat lon" ;
##     float time()    ; 
##         time:standard_name = "time" ;
##         time:long_name = "time of measurement" ;
##         time:units = "days since 2000-01-01 00:00:00" ;
##     char gages(linkDim,IDLength)    ; 
##         gages:units = "usgsId" ;
##         gages:_FillValue = "               " ;
## 
## // global attributes (2):
##     :featureType = "point"
##     :history = "Created Thu Sep 03 09:30:05 2015"
print(ncdump(newFile, 'gages', quiet=TRUE))
##   [1] "           g255" "           g266" "           g213"
##   [4] "           g253" "           g277" "           g208"
##   [7] "           g259" "           g240" "           g241"
##  [10] "           g257" "           g217" "           g216"
##  [13] "           g238" "           g228" "           g247"
##  [16] "           g229" "           g231" "           g258"
##  [19] "           g069" "           g230" "           g061"
##  [22] "           g244" "           g063" "           g222"
##  [25] "           g263" "           g262" "           g066"
##  [28] "           g256" "           g252" "           g115"
##  [31] "           g141" "           g140" "           g225"
##  [34] "           g103" "           g067" "           g065"
##  [37] "           g091" "           g090" "           g260"
##  [40] "           g187" "           g227" "           g146"
##  [43] "           g112" "           g210" "           g264"
##  [46] "           g133" "           g132" "           g163"
##  [49] "           g226" "           g128" "           g221"
##  [52] "           g092" "           g215" "           g110"
##  [55] "           g177" "           g142" "           g139"
##  [58] "           g234" "           g199" "           g164"
##  [61] "           g106" "           g123" "           g169"
##  [64] "           g118" "           g209" "           g170"
##  [67] "           g105" "           g200" "           g197"
##  [70] "           g126" "           g223" "           g186"
##  [73] "           g153" "           g120" "           g147"
##  [76] "           g175" "           g107" "           g179"
##  [79] "           g224" "           g096" "           g124"
##  [82] "           g154" "           g089" "           g151"
##  [85] "           g211" "           g176" "           g203"
##  [88] "           g138" "           g167" "           g099"
##  [91] "           g125" "           g155" "           g152"
##  [94] "           g150" "           g181" "           g207"
##  [97] "           g161" "           g097" "           g190"
## [100] "           g058" "           g280" "           g278"
## [103] "           g050" "           g144" "           g174"
## [106] "           g098" "           g188" "           g059"
## [109] "           g282" "           g281" "           g055"
## [112] "           g113" "           g204" "           g104"
## [115] "           g166" "           g165" "           g189"
## [118] "           g121" "           g056" "           g052"
## [121] "           g178" "           g172" "           g171"
## [124] "           g300" "           g295" "           g101"
## [127] "           g191" "           g094" "           g093"
## [130] "           g316" "           g122" "           g057"
## [133] "           g180" "           g049" "           g303"
## [136] "           g302" "           g325" "           g095"
## [139] "           g317" "           g314" "           g343"
## [142] "           g053" "           g212" "           g040"
## [145] "           g341" "           g100" "           g194"
## [148] "           g087" "           g082" "           g332"
## [151] "           g330" "           g198" "           g195"
## [154] "           g319" "           g086" "           g083"
## [157] "           g235" "           g201" "           g326"
## [160] "           g318" "           g085" "           g205"
## [163] "           g108" "           g196" "           g320"
## [166] "           g088" "           g310" "           g116"
## [169] "           g145" "           g304" "           g109"
## [172] "           g074" "           g232" "           g068"
## [175] "           g192" "           g184" "           g183"
## [178] "           g084" "           g305" "           g296"
## [181] "           g071" "           g102" "           g070"
## [184] "           g292" "           g288" "           g315"
## [187] "           g313" "           g290" "           g182"
## [190] "           g080" "           g078" "           g173"
## [193] "           g299" "           g297" "           g168"
## [196] "           g162" "           g289" "           g157"
## [199] "           g156" "           g079" "           g309"
## [202] "           g308" "           g143" "           g075"
## [205] "           g298" "           g136" "           g294"
## [208] "           g285" "           g248" "           g119"
## [211] "           g117" "           g076" "           g137"
## [214] "           g261" "           g291" "           g322"
## [217] "           g321" "           g287" "           g283"
## [220] "           g249" "           g329" "           g242"
## [223] "           g336" "           g111" "           g073"
## [226] "           g038" "           g324" "           g323"
## [229] "           g032" "           g286" "           g251"
## [232] "           g345" "           g254" "           g340"
## [235] "           g051" "           g338" "           g337"
## [238] "           g062" "           g339" "           g328"
## [241] "           g327" "           g293" "           g352"
## [244] "           g348" "           g347" "           g346"
## [247] "           g361" "           g054" "           g021"
## [250] "           g243" "           g018" "           g047"
## [253] "           g334" "           g333" "           g265"
## [256] "           g335" "           g359" "           g006"
## [259] "           g356" "           g355" "           g033"
## [262] "           g127" "           g378" "           g376"
## [265] "           g369" "           g239" "           g366"
## [268] "           g365" "           g362" "           g360"
## [271] "           g358" "           g357" "           g349"
## [274] "           g214" "           g077" "           g012"
## [277] "           g311" "           g048" "           g301"
## [280] "           g039" "           g034" "           g030"
## [283] "           g025" "           g276" "           g019"
## [286] "           g273" "           g268" "           g267"
## [289] "           g131" "           g379" "           g377"
## [292] "           g375" "           g374" "           g245"
## [295] "           g371" "           g114" "           g368"
## [298] "           g367" "           g237" "           g364"
## [301] "           g363" "           g354" "           g353"
## [304] "           g350" "           g284" "           g014"
## [307] "           g081" "           g193" "           g031"
## [310] "           g312" "           g046" "           g041"
## [313] "           g036" "           g035" "           g159"
## [316] "           g028" "           g027" "           g023"
## [319] "           g022" "           g020" "           g275"
## [322] "           g015" "           g270" "           g269"
## [325] "           g388" "           g382" "           g381"
## [328] "           g380" "           g373" "           g372"
## [331] "           g236" "           g233" "           g016"
## [334] "           g351" "           g220" "           g218"
## [337] "           g342" "           g386" "           g064"
## [340] "           g060" "           g307" "           g306"
## [343] "           g045" "           g042" "           g037"
## [346] "           g160" "           g389" "           g029"
## [349] "           g026" "           g149" "           g274"
## [352] "           g017" "           g272" "           g013"
## [355] "           g008" "           g135" "           g005"
## [358] "           g387" "           g130" "           g385"
## [361] "           g384" "           g383" "           g250"
## [364] "           g246" "           g370" "           g391"
## [367] "           g219" "           g344" "           g399"
## [370] "           g206" "           g331" "           g202"
## [373] "           g072" "           g395" "           g394"
## [376] "           g185" "           g393" "           g044"
## [379] "           g043" "           g129" "           g390"
## [382] "           g158" "           g024" "           g279"
## [385] "           g148" "           g403" "           g402"
## [388] "           g401" "           g400" "           g271"
## [391] "           g398" "           g397" "           g396"
## [394] "           g011" "           g010" "           g009"
## [397] "           g392" "           g007" "           g134"
## [400] "           g004" "           g003" "           g002"
## [403] "           g001"

Run the model with doubled precip

Set the symlink to the new Route_Link.nc file.

jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek/DOMAIN> l
lrwxrwxrwx  2 jamesmcc rap   93 Sep  9 09:42 Route_Link.nc -> 
              ../../DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc

Set up several model parameters in the namelists.

namelist.hrldas:

 START_YEAR  = 2012
 START_MONTH  = 08
 START_DAY  = 01
 START_HOUR  = 00
 START_MIN   = 00
 KHOUR = 48
 FORCING_TIMESTEP = 3600
 NOAH_TIMESTEP    = 900
 OUTPUT_TIMESTEP  = 7200
 FORC_TYP = 4

hydro.namelist:

 out_dt = 15 ! minutes
 DTRT_CH = 60
 DTRT_TER = 6
 SUBRTSWCRT = 1
 OVRTSWCRT = 1
 CHANRTSWCRT = 1
 channel_option = 2
 route_link_f = "DOMAIN/Route_Link.nc"

Run without nudging and using precip doubling. The wrf_hydro.noNudging_doublePrecip.exe binary was built with the following environment variables: PRECIP_DOUBLE=1 WRF_HYDRO_NUDGING=0 using the daBranch of wrf_hydro_model.

jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.noNudging_doublePrecip.exe
. . .

jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmo
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmo

Transform output to timeSlices

Check chanobs and frxst pts dimensions

jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek/run.pmo> wc -l frxst_pts_out.txt
77779 frxst_pts_out.txt  ## = 403stns*(48hours*4output/hr+1output@0)
frxst <- ReadFrxstPts("~/WRF_Hydro/Col_Bldr_Creek/run.pmo/frxst_pts_out.txt", stId='character')
ggplot2::ggplot(subset(frxst, POSIXct < as.POSIXct('2012-08-01 06:00:00', tz='UTC')), ## all the action is early on
                ggplot2::aes(x=POSIXct, y=q_cms, color=st_id)) +
  ggplot2::geom_line()  +
  ggplot2::scale_color_discrete(guide='none') +
  ggplot2::theme_bw()

Plot it by order

rlOrder <- ncdump(rlFile, 'order', quiet=TRUE)
names(rlOrder) <- paste0('g',formatC(rlLink,width=3,flag='0'))
frxst$order <- rlOrder[trimws(frxst$st_id)]
ggplot2::ggplot(subset(frxst, POSIXct < as.POSIXct('2012-08-01 06:00:00', tz='UTC')), 
                ggplot2::aes(x=POSIXct, y=q_cms, color=st_id)) +
  ggplot2::geom_line()  +
  ggplot2::scale_color_discrete(guide='none') +
  ggplot2::facet_wrap(~order, ncol=1, scales='free_y') +
  ggplot2::theme_bw()

Create timeslice files

chanObsFiles <- list.files('~/WRF_Hydro/Col_Bldr_Creek/run.pmo/', pattern = 'CHANOBS_DOMAIN', full.names=TRUE)
sliceResolutionMin <- 15
outputDir <- '~/WRF_Hydro/Col_Bldr_Creek/run.pmo/timeSlicePMO/'
dir.create(outputDir)
sliceFiles <- ChanObsToTimeSlice(chanObsFiles, sliceResolutionMin, outputDir) 

Put these in place

jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> rm nudgingTimeSliceObs 
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ln -sf run.pmo/timeSlicePMO/ nudgingTimeSliceObs

Create the nudging parameter file

min(ncdump(rlFile, 'Length', quiet=TRUE))  ## use R less than this /2.
## [1] 50
gageParams <- data.frame( gageId = paste0('g',formatC(1:403,width=3,flag='0')), stringsAsFactors=FALSE)
gageParams$R=24
gageParams$G=1
gageParams$tau=20
MkNudgingParams(gageId=gageParams$gageId, R=gageParams$R, 
               G=gageParams$G, tau=gageParams$tau, 
               outFile='~/WRF_Hydro/Col_Bldr_Creek/run.pmo/nudgingParams.PMOallGages.nc', 
               overwrite=TRUE)
## File: ~/WRF_Hydro/Col_Bldr_Creek/run.pmo/nudgingParams.PMOallGages.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (2):
##     stationIdStrLen = 15 ; 
##     stationIdInd = UNLIMITED ; // (403 currently)
## variables (4):
##     char stationId(stationIdInd,stationIdStrLen)   (Chunking: [15,1])   ; 
##         stationId:units = "-" ;
##         stationId:long_name = "USGS station identifer" ;
##     float R(stationIdInd)   (Chunking: [1])   ; 
##         R:units = "meters" ;
##         R:long_name = "Radius of influence in meters" ;
##     float G(stationIdInd)   (Chunking: [1])   ; 
##         G:units = "-" ;
##         G:long_name = "Amplitude of nudging" ;
##     float tau(stationIdInd)   (Chunking: [1])   ; 
##         tau:units = "minutes" ;
##         tau:long_name = "Time tapering parameter half window size in minutes" ;

Then put the param file in place.

jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ln -s run.pmo/nudgingParams.PMOallGages.nc nudgingParams.nc

Assimilate the observations using regular forcing

The wrf_hydro.nudging.exe binary was built with the following environment variables: PRECIP_DOUBLE=0 and WRF_HYDRO_NUDGING=1 using the daBranch of wrf_hydro_model.

jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.nudging.exe
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmoNudge
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmoNudge

Run the model without nudging using regular forcingings.

The wrf_hydro.noNudging.exe binary was built with the following environment variables: PRECIP_DOUBLE=0 and WRF_HYDRO_NUDGING=0 using the daBranch of wrf_hydro_model.

jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.noNudging.exe
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmoNoNudge
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmoNoNudge

Visualize

77779 frxst_pts_out.txt ## = 403stns*(48hours*4output/hr+1output@0)

frxstNoNudge <- ReadFrxstPts("~/WRF_Hydro/Col_Bldr_Creek/run.pmoNoNudge/frxst_pts_out.txt", stId='character')
frxstNudge <- ReadFrxstPts("~/WRF_Hydro/Col_Bldr_Creek/run.pmoNudge/frxst_pts_out.txt", stId='character')

## make sure the abscissae columns are equal before taking the differences
if(all(plyr::laply(NamedList(names(frxst)), function(colName) all(frxst[[colName]]==frxstNoNudge[[colName]]))[-6:-7])) 
{
  forcErr <- frxst
  forcErr$`q err obs-model (cms)` <- frxst$q_cms - frxstNoNudge$q_cms
  forcErr$`q err obs-model (cfs)` <- frxst$q_cfs - frxstNoNudge$q_cfs
  forcErr$q_cms <- forcErr$q_cfs <- NULL
  forcErr$error <- 'forcing'
} else warning('Abscissae do not line up with frxstNoNudge, please investigate', immediate.=TRUE)

if(all(plyr::laply(NamedList(names(frxst)), function(colName) all(frxst[[colName]]==frxstNudge[[colName]]))[-6:-7])) 
{
  nudgeErr <- frxst
  nudgeErr$`q err obs-model (cms)` <- frxst$q_cms - frxstNudge$q_cms
  nudgeErr$`q err obs-model (cfs)` <- frxst$q_cfs - frxstNudge$q_cfs
  nudgeErr$q_cms <- nudgeErr$q_cfs <- NULL
  nudgeErr$error <- 'nudging'
} else warning('Abscissae do not line up with frxstNudge, please investigate', immediate.=TRUE)

err <- rbind(forcErr, nudgeErr)

ggplot2::ggplot(subset(err, POSIXct < as.POSIXct('2012-08-01 06:00:00', tz='UTC')), ## all the action is early on
                ggplot2::aes(x=POSIXct, y=`q err obs-model (cms)`, color=st_id)) +
  ggplot2::geom_line()  +
  ggplot2::scale_color_discrete(guide='none') +
  ggplot2::facet_wrap(~error, ncol=1) + 
  ggplot2::theme_bw()